ANALYSIS OF FULLY DISCRETE FINITE ELEMENT METHODS FOR A SYSTEM 
OF DIFFERENTIAL EQUATIONS MODELING SWELLING DYNAMICS OF 

POLYMER GELS 

XIAOBING FENG* AND YINNIAN HEt 

Abstract. The primary goal of this paper is to develop and analyze some fully discrete finite element methods for 
a displacement-pressure model modeling swelling dynamics of polymer gels under mechanical constraints. In the model, 
the swelling dynamics is governed by the solvent permeation and the elastic interaction; the permeation is described by a 
pressure equation for the solvent, and the elastic interaction is described by displacement equations for the solid network 
of the gel. The elasticity is of long range nature and gives effects for the solvent diffusion. It is the fluid-solid interaction in 
the gel network drives the system and makes the problem interesting and difficult. By introducing an "elastic pressure" (or 
"volume change function" ) we first present a reformulation of the original model, we then propose a time-stepping scheme 
which decouples the PDE system at each time step into two sub-problems, one of which is a generalized Stokes problem for 
the displacement vector field (of the solid network of the gel) and another is a diffusion problem for a "pseudo-pressure" field 
(of the solvent of the gel). To make such a multiphysical approach feasible, it is vital to find admissible constraints to resolve 
the uniqueness issue for the generalized Stokes problem and to construct a "good" boundary condition for the diffusion 
equation so that it also becomes uniquely solvable. The key to the first difficulty is to discover certain conservation laws 
(or conserved quantities) for the PDE solution of the original model, and the solution to the second difficulty is to use the 
generalized Stokes problem to generate a boundary condition for the diffusion problem. This then lays down the theoretical 
foundation for one to utilize any convergent Stokes solver (and its code) together with any convergent diffusion equation 
solver (and its code) to solve the polymer gel model. In the paper, the Taylor-Hood mixed finite element method combined 
with the continuous linear finite element method are chosen as an example to present the ideas and to demonstrate the 
viability of the proposed multiphysical approach. It is proved that, under a mesh constraint, both the proposed semi-discrete 
(in space) and fully discrete methods enjoy some discrete energy laws which mimic the differential energy law satisfied by 
the PDE solution. Optimal order error estimates in various norms are established for the numerical solutions of both the 
semi-discrete and fully discrete methods. Numerical experiments are also presented to show the efficiency of the proposed 
approach and methods. 

Key words. Gels, soft matters, poroelasticity. Stokes equations, finite element methods, inf-sup condition, fully discrete 
schemes, error estimates. 
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1. Introduction. A gel is a soft poroelastic material which consists of a solid network 
and a colloidal solvent. The solid network spans the volume of the solvent medium. The 
solvent can permeate through the solid network and the permeation can be controlled by 
external forces. Both by weight and volume, gels are mostly liquid in composition and 
thus exhibit densities similar to liquids. However, they have the structural coherence of a 
solid and can be deformed. A gel network can be composed of a wide variety of materials, 
including particles, polymers and proteins, which then gives different types gels such 
hydrogels, organogels and xerogels (cf. Gels have some fascinating properties, in 

particular, they display thixotropy which means that they become fluid when agitated, 
but resolidify when resting. In general, gels are apparently solid, jelly-like materials, they 
exhibit an important state of matter found in a wide variety of biomedical and chemical 
systems (cf. |9l UHl UHl [20] and the references therein). 

This paper develops and analyzes some fully discrete finite element methods for a 
displacement-pressure model for polymer gels. The model, which was proposed by M. 
Doi et al in [HI [IHl 120], describes swelling dynamics of polymer gels (under mechanical 
constraints). Let C M'' (o? = 1, 2, 3) be a bounded domain and denote the initial region 
occupied by the gel. Let u(a;, t) denote the displacement of the gel at the point x & VL 
in the space and at the time t, Vs(x,t) and p{x,t) be the velocity and the pressure of 
the solvent at {x,t). Following |9], the governing equation for the swelling dynamics of 
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polymer gels are given by 



(1.1) 
(1.2) 
(1.3) 



div (cr(u) — pi) 
^(v, - Ut) 
div {(f)Ut + (1 - 0)v,) 



0)Vp, 



0. 



Here ^ is the friction constant associated with the motion of the polymer relative to the 
solvent, is the volume fraction of the polymer, I denotes the d x d identity matrix, 
and cr(u) stands for the stress tensor of the gel network, which is given by a constitutive 
equation. In this paper, we use the following linearized form of the stress tensor: 



1.4) 



{K- -G')divu/ + 2Ge(u) 
3 



1 



£(u):=-(Vu + Vu^), 



where K and G are respectively the bulk and shear modulus of the gel (cf. [HI E]). We 
remark that (1.1 ) stands for the force balance, ( |1.2 ) states Darcy's law for the permeation 
of solvent through the gel network, and (1.3) describes the incompressibility condition. 
In addition, if we introduce the total stress a{u,p) := cr(u) — pi, then equation (1.1) 
becomes diva(u,p) = 0. 



Substituting (1.4) into (1.1) and (1.2) into (1.3) yield the following basic equations for 



swelling dynamics of polymer gels (see [19] ) 



(1.5) 
(1.6) 



aVdivu + (3Au = Vp, 
divut = nAp, 



a:=K 



G 



/3:=G, 



which hold in the space-time domain Qt '■= ^ x (0,T) for some given T > 0. 

To close the above system, we need to prescribe boundary and initial conditions. Only 
one initial condition is required for the system, which is 



;i-7) 



u(-,0) = uo( 



in Q. 



Various sets of boundary conditions are possible and each of them describes a certain type 
mechanical condition and solvent permeation condition (cf. [T9l 120]). In this paper we 
consider the following set of boundary conditions 



(1.^ 



(cr(u) -pl)u = f. 



dp 
du 



onriT ■.= dnx {0,T), 



where i' denotes the outward normal to dfl 
applied on the boundary of the gel. Since 



(01 



means that the mechanical force f is 
, hence, (1.8)2 implies that 



- Ut 



the solvent can not permeate through the gel boundary. We also remark that the force 
function f must satisfy the compatibility condition 



fdS = 0. 



an 



Problem (1.5)-(1.8) is interesting and difficult due to its multiphysical nature which 



describes the complicate fluid and solid interaction inside the gel network. It is numerically 
tricky to solve because it is difficult to design a good and workable time-stepping scheme. 
For example, one natural attempt would be at each time step first to solve a Poisson 
problem for p and then to solve a linear elasticity problem for u. However, this strategy 
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is difficult to realize because there is no good way to compute the source term div ut for 
the Poisson equation. In fact, the strategy even has a difficulty to start due to the fact 
that no initial condition is provided for the pressure p. 



To overcome the difficulty, in this paper we shall use a reformulation of system (1.5)- 



(1.6), which is now introduced. Define 

(1.9) q := divu. 

Physically, q measures the volume change of the solid network of the gel, and often called 



■'elastic pressure" or "volume change function". Taking divergence on (1.5) yields 

{a + (3)Aq = Ap, 



which and (1.6) imply that q satisfies the following diffusion equation 



(1.10) qt = DAq, D:=k{K+^G). 

However, the usefulness of the above diffusion equation is hampered by the lack of bound- 
ary condition for q. We like to note that the above diffusion equation for q was ffist noticed 
by M. Doi [S], but it was not utilized before exactly because of the lack of boundary con- 



dition for q. Nevertheless, using the new variable q we can rewrite (1.5)-(1.6) as 



^1.11) /3Au = Vp, p := p — aq, 

1.12) divu = g, 

1.13) qt = nAp, p = p + aq. 



An immediate consequence of the above reformulation is that (1.11)-(1.12) implies {u,p) 
satisfies the generalized Stokes equations with q being the source term at each time t, and 
q satisfies a diffusion equation and it interacts with (u,p) only at the boundary dQ. 

This is a key observation because it not only reveals the underlying physical process 
of swelling dynamics of the gel, but also gives the "right" hint on how the problem should 
be solved numerically. This indeed motivates the main idea of this paper, that is, at 
each time step, we first solve the generalized Stokes problem for {u,p), which in turn 
provides (implicitly) an updated boundary condition for q , we then use this new boundary 
condition to solve the diffusion equation for q. The process is repeated iteratively until 
the final time step is reached. However, in order to make this idea work, there is one 
crucial issue needs to be addressed. That is, for a given q the generalized Stokes problem 
for (u,p) is only unique up to additive constants. Clearly, how to correctly enforce the 
uniqueness of the generalized Stokes problem is the bottleneck of this approach. It is 
easy to understand that one can not use arbitrary constraints to fix (u,p) because this 
will lead to bad or even divergent numerical schemes if the exact PDE solution does not 
satisfy the constraints. Instead, the constraints which can be used to fix (u, p) should be 
those satisfied by the exact solution of the PDE system. To the end, we need to discover 
some invariant (or conserved) quantities for the exact PDE solution. It turns out that 
the situation is precisely what we anticipated and wanted. We are able to show that the 
exact PDE solution (u, p, p, q) satisfies the following identities (see Section |2] below for a 
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proof): 








(1.14) 


/ Qix 
1 

Jn 


t)dx 


y 


(1.15) 


/ p{x 
Jn 


t)dx 




(1.16) 


/ u(x,t)- 
Jan 


V dS 





n 

CdCn 



qo{x)dx := / divuo{x)dx, 



n 



1 
d 

uo(x) 



f{x, t) ■ X dS, 



where d denote the dimension of Q and 



(1.17) 



Crf := a + 



d 



on 



-fv -r g 



an 

V dS 



iid 
iid 



C, 



11 



2, 
3. 



Obviously, the right-hand sides of (1.14) and (1.16) are constants. The right-hand side 



of (1.15) is also a constant provided that f is independent of t, otherwise, it is a known 



function of t. In this paper, we shall only consider the case that f is independent of t. It 



follows from (|1.14|) and (|1.15|) that 
(1.18) 



I p{x, t)dx = Cp := Cp - aCg = - k [ ■ 



f(x) ■ X dS. 



It is clear now that (1.18) and (1.16) provide two natural conditions which can be used 



to uniquely determine the solution (u,p) to the generalized Stokes problem ( 1.11 )-( 1.12) 



for a given source term q. This then leads to the following time-discretization for problem 



(1.5)-(1.8): 



Algorithm 1: 

(i) Set g° = go '■= divuo and u° := uq. 

(ii) For n = 0, 1, 2, ■ ■ ■ , do the following two steps 
Step 1: Solve for (u""''^,^"'"'"^) such that 



(1.19) 




= 0, 


in Qx, 


(1.20) 


divu"+i 


= g" 


in Qt, 


(1.21) 


(3— 


= f 


on dri 


(1.22) 


(p"+M) = C^, (u"+\z/) 







Step 2: Solve for such that 



(1.23) 

(1.24) 
(1.25) 

where dtq"'^^ 



0, 



a 



dq 



n+l 



rn+1 



dv 

(g"+\ 1) 



dv 



in 

on OVLt) 



C, 



We note that (1.23) is the implicit Euler scheme, which is chosen just for the ease of pre- 



sentation, it can be replaced by other time-stepping schemes. (1.24) provides a Neumann 
boundary condition for g"+^. Another subtle issue is the role which the initial value uq 
plays in the algorithm. Seems Uq is only needed to produce go and there is no need to have 
u° in order to execute the algorithm. However, to ensure the stability and convergence 
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of the algorithm, it turns out that u° not only needs to be provided but also must be 
carefully constructed when the algorithm is discretized (see Sections [3] and |4]). 

The above algorithm has a couple attractive features. First, it is easy to use. Second, it 
allows one to make use of any available numerical methods (finite element, finite difference, 
finite volume, spectral and discontinuous Galerkin) and computer codes for the Stokes 
problem and the Poisson problem to solve the gel swelling dynamics model (1.5)-(1.8). We 
remark that an almost same model as ( 1.5[ )-( |L8 ) also arise from different applications in 
poroelasticity and soil mechanics and are known as Bolt's consolidation model (cf. [21 IH] 
and the references therein). [2] proposed and analyzed a standard finite element method 
which directly approximates (u, p) under the (restrictive) divergence-free assumption on 
the initial condition uq. 

This paper consists of four additional sections. In Section[2| we first introduce notation 
used in this paper. We then present a PDE analysis for the gel swelling dynamics model 
(1.5)-(1.8), which includes deriving a dissipative energy law, establishing existence and 



uniqueness, and in particular, proving the conservation laws stated in ( 1.14 )-( 1.16) and 
(1.18). In Section [3| we first propose a semi-discrete (in space) finite element discretization 
for problem (1.5)-(1.8) based on the multiphysical reformulation ( 1.11 )-( 1.13). The well- 



known Taylor-Hood mixed element and the Pi conforming finite element are used as an 
example to present the ideas. It is proved that the solution of the semi-discrete method 
satisfies a discrete energy law which mimics the differential energy law enjoyed by the 
PDE solution, and the semi-discrete numerical solution also satisfies the conservation 



laws (1.14)-(1.16) and (1.18). We then derive optimal order error estimates in various 



norms for the semi-discrete numerical solution. In Section |4| fully discrete finite element 
methods are constructed by combining the time-stepping scheme of Algorithm 1 and the 
semi-discrete finite element methods of Section IH The main results of this section include 
proving a fully discrete energy law for the numerical solution and establishing optimal 
order error estimates for the fully discrete finite element methods. Finally, in Section [5| 
we present some numerical experiments to gauge the efficiency of the proposed approach 
and methods. 



2. PDE analysis for problem ( |1.5[ )-pr8|. The standard Sobolev space notation is used 
in this paper, we refer to [H Ej [18j for their precise definitions. In particular, (-, ■) and 
(■, ■) denote respectively the standard L'^{Q) and L'^{dfl) inner products. For any Banach 
space B, we let B = [BY and use B' to denote its dual space. In particular, we use (■, ■)duai 
and (-, ■)duai to denote the dual products on {H^{Q)y x H^{Q) and H^2[dQ) x H2[dfl), 
respectively. || • \\lp{b) is a shorthand notation for || ■ ||lp({o,t);B)- 

We also introduce the function spaces 

LliQ) := {q G L^Q); (g, 1) = 0}, X := {v G H^Q); (v, u) = 0}. 

It is well known [18J that the following so-called inf-sup condition holds in the space 
X X Ll{n): 



(2.1) 



sup 
vex 



(divv,v?) 



Vv 



>ao\\f\\L^ \f!feLl{Q), ao>0. 



Throughout the paper, we assume f2 C M'^ be a bounded polygonal domain such that 
A : Hq{Q) n H'^{Q) Li^iyt) is an isomorphism; see [TTl [12]. In addition, C is used 
to denote a generic positive constant which is independent of u, p, p, g, and the mesh 
parameters h and At. 



We now give a definition of weak solutions to (1.5)-(1.8). 
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Definition 2.1. Let (uo,f) e U\n) x U'^dn), and (f , l)duai = 0. Given T >0, a 
tuple {u,p) with 



is called a weak solution to (1.5)-(1.8), if there hold for almost every t G [0,T] 

(2.2) ((divu)i,<^)^^^j + /c(Vp,V<^) =0 \/veH\Q), 

(2.3) a(divu,divv) +/5(Vu,Vv) - (p,divv) = (f,v)duai Vv G H^(fi), 

(2.4) u(0) = uo. 



Similarly, we define weak solutions to problem (1.11)-(1.13), (1.7)-(1.8). 



Definition 2.2. Let (uo,f) G x and (f, l)duai = 0. Given T >0, a 

triple (u, p, q) with 

u G L°°(0,T;Hi(r])), pG L\Q,T-L\n)), 

q G L°°(0,^;L2(^])) nif^(0,r;i/-i(fi)), aq + pe L\0,T; H\n)), 



is called a weak solution to (1.11 )-(1.13), (1.7)-(1.8) if there hold for almost every t & [0;^] 
(2.5) /?(Vu, Vv) - (p,divv) = (f,v)duai Vv G Hi(fi), 



(2.6) 
(2.7) 
(2.8) 
(2.9) 



(divu,v?) = {q,(p) 
u(0) = Uo, 

(^*'^)dual + ^(^("^ + ^'^^) =0 

g(0) = go := divuo. 



G L2(fi), 
G i/^(fi). 



Remark 2.1. (^aj Glearly, p := aq + p gives back the pressure in the original formula- 
tion. What interesting is that both q andp are only L'^ -functions in the spatial variable but 
their combination aq + p is an -function. In other words, the new formulation provides 
an — decomposition for the pressure p. It turns out that this decomposition will have 
a significant numerical impact because it allows one to use low order (hence cheap) finite 
elements to approximate q and p but still to be able to approximate the pressure p with 
high accuracy. 



(b) (2.8) implicitly imposes the following boundary condition for q: 



(2.10) 



a 



dq 
du 



dp 
dv 



on dQ. 



Since problem (2.2)-(2.4) consists of two linear equations, its solvability should fol- 



lows easily if we can establish a priori energy estimates for its solutions. The following 
dissipative energy law just serves that purpose. 



Lemma 2.3. Every weak solution (u,p) of problem (1.5)-(1.8) satisfies the following 
energy law: 



(2.11) 



E{t) + K \\ Vp{s 



Il2 ds 



E{0) VtG(0,T], 
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where 

Moreover, 
(2.12) 



E{t) := ^ [/?|| Vu(t) 11^2 + all divu(t) - 2(f , u(t))duai 



(divu)t ||l2(h-i) = «;|| ||l2(h-i) < -E(0)2. 



Proof. We first consider the case u G H^{0,T;li^{Q)). Setting (f = p in (2.2) and 
V = in (2.3) yield 

[div ut{t), p{t)) + 4 Vp\\l, = 0, 



d 

dt 



flldivu(t) ||2, + ^||Vu(t) Hi. 



d 
dt 



(f,u(t))duai+ (p(^),divut(t)). 



(2.11 ) follows from adding the above two equations and integrating the sum in t over the 



interval (0, s) for any s G (0,T]. 



If (u,p) is only a weak solution, then could not be used as a test function in (2.3). 



However, this difficulty can be easily overcome by using as a test function in (2.3) 



where u'' denotes a mollification of u through a symmetric mollifier (cf. [HI Chapter 7]) 
and by passing to the limit 6^0. 



Finally, (2.12) follows immediately from (2.2) and (2.11). The proof is completed. □ 



(1.11 


)-( 


1.13 


), ( 


1.7 


)-( 


1.8 



now is replaced by the following equivalent energy law: 

(2.13) J{t) + K [ \\V[p{s) + aq{s)]\\l2ds = J{0) VtG(0,T], 

Jo 



and (2.12) is replaced by 

II qt Wl^h-^) < v^ll V{p + aq) ||l2(l2) < J(0)i 



(2.14) 

Where 



J{t) 



Vu(t) 11^2 + all g(t) Ilia - 2(f , u(t))dual 



We also note that the weak solution to problem (1.11 )-( fl.l3 ), ( 1.7)-( |L8| is understood 



the sense of Definition 2^ 



Lemma 2.4. Weak solutions of proble m (|L5| -( pr8|) and problem ( |1.11[ )-( [lI3| ),([L7)) 



(|1.8) satisfy the conservation laws (1.14)-(1.16), and (1.18). 

Proof. (1.14) and (1.16) follows immediately from taking if) = 1 in ( |2.8[ ), = 1 in 



both (2.6) and (2.2) followed by integrating in t and appealing to the divergence theorem. 



To prove (1.15), taking v = x in (2.3) and using the identities Vx = I and divx = d 
we get 



a 



(divu,rf) +/3(Vu,/) - [p,d) = (f,x)duai 



Hence, 



{p, l) = a (div u, l) + ^ (div u, l) - ^ (f , x)duai 
a + ^) (divuo, l) - ^(f,x)duai 



CdCq — ^(f,a^)duah 
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where q and Cq are defined in (1.17) and (1.14). So we obtain (1.15). 

Similarly, (1.18) follows by taking v = x in (2.5). The proof is complete. □ 



With the help of the above two lemmas, we can show the solvability of problem (1.5)- 



Theorem 2.5. Suppose (uo,f) G H^(fi) x H~^(5f2), and (f, l)duai = 0. Then there 
exists a unique weak solution to ( 1.5[ )-(1.8) in the sense of Definition 2.1, and there exists 
a unique solution to to ( |1.11[ )-(1.13), (1.7)-(1.8). 

Proof. The uniqueness is an immediate consequence of the energy laws (2.11) and 



(2.13) and the conservation laws (1.14)-(1.16) and (1.18). 

The existence can be easily proved by using (abstract) Galerkin method and the 
standard compactness argument (cf. [IB]). The energy laws (2.11) and (2.13) provide the 
necessary uniform estimates for the Galerkin approximate solutions. We omit the details 
because the derivation is standard. □ 

Again, exploiting the linearity of the system, we have the following regularity results 
for the weak solution. 

Theorem 2.6. Let {u,p,p,q) be the weak solution of (2.5)-(2.9) Then there holds 
the following estimates: 

(2.15) 
(2.16) 
(2.17) 
(2.18) 
(2.19) 
(2.20) 



(2.21) 



ViVut ||l2(^2) + ^/a\\ Viqt ||l2(l2) + v^|| VtVp ||Loo(i2) < J(0)2. 
tVut ||l=°(l2) + Va\\ tqt + v^|| tVpt < [2J(0)]i 

WtquhHH-^) < [2J(0)]i 
K^/a\\ A/tAp 11^2(^2) = ^/a\\ Viqt < -^(0)i 

/tv/a||tAp 11^00(^2) = i|L-(L2) < [2J(0)]i 

II ||l2(l2) + II Vp^||Loo(i2) + a^/K\\ Vg ||l2(l2) 

+«v^|| ViVq |Uoo(i2) < C{Cp, J(0)). 

P\\ AU |U^(^2) = II Vp |Uoo(i2) < C{Cp, J(0)). 



Where C(T,Cp, J(0)) denotes some positive constant which depends on T,Cp and J(0). 

Proof. Differentiating each of equations (2.5), (2.6) and (2.8) with respect to t yields 
(note that f is assumed to be independent of t) 

(2.22) /3(Vuf, Vv) - {pt, div v) =0 Vv G H^(fi), 

(2.23) (divut,<^) = {quip) V<^ G L\n), 

(2.24) (g«,^)d,,i + '^(V(ag + p)t,VV^) =0 W^p e H\n). 



Taking v = tut in (2.22), ip = tpt in (2.23), = tpt = t{p + aq)t in (2.8), and adding the 
resulted equations we get 



Vui 11^2+ toll 11^2 + ^^(t||Vp 11^2 



K , 



Vp 



L2- 



Integrating in t over (0, T) gives (2.15). 

Alternatively, setting v = t'^Uu in (2.22), ip = t^pt in (2.23) (after differentiating in t 
one more time), ip = t^pt = t^{j) + aq)t in (2.24), and adding the resulted equations we 
have 



1 d 
2dt 



t^PW Vut + t^a\\ qt 11^2 + t^4 Vpt 11^2 



Vut||i2 



^11 (It 



lL2- 
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Integrating in t over (0,T) and using (2.15) yield (2.16). 



(2.17)-(2.19) are immediate consequences of (2.24), (2.15), and (2.16). 



Choose 00 £ C^{fi) such that (0o, 1) = Cp- Then p 
version of the inf-sup condition (cf. f^) we have 



G Lq{Q). By an equivalent 



V(p - 0o) IU2 < sup 
veHi(n) 



(V(p-0o),v) 



Vv 



L2 



^ ^ -/5(Vu,Vv) 

< q;q sup - 

vGHi{Q) 



< a. 



-1 



I Vv 11^2 
Vu 11^2 + v^ll 00 IU2 



ttQ ^v^ll 01 







Hence, 



Vp ||l2 < II V0O ||l2 + a~Q 



H Vu||i2 + v^||0o|U2 
which together with (2.13) and (2.15) imply that (recall p = p + aq) 

II Vp 11^2(^2) < v^ll Vp |Uoc(i2) < Vf\\ V0O IU2 + a^^Vfbj{Q)-2 + v^ll 00 IU2 



oi\f^\\ Vg ||L2(i2) < J(0)^ + v^ll V0O ||l2 + a^^V^ /5J(0)^ + v^|| 0o IU2 



av^ll ViWq |Uoc(i2) < J(0)^ + v^ll V0O ||l2 + tto /?^(0)^ + v^|| 0o lU^ 



Hence, i pM ) holds. 

Finally, (2.21) follows immediately from the equation f3Au = Vp and (2.20). The 
proof is complete. □ 

3. Semi-discrete finite element methods in space. The goal of this section is to present 
the ideas and specific semi-discrete finite element methods for discretizing the variational 
problem (2.5)-(2.9) in space based on the multiphysical (deformation and diffusion) ap- 
proach. That is, we shall approximate (u,p) using a (stable) Stokes solver and approx- 
imate g by a convergent diffusion equation solver. The combination of the Taylor-Hood 
mixed finite element [1] and the conforming Pi finite element is chosen as a specific ex- 
ample to present the ideas. 

3.1. Formulation of finite element methods. Assume Q G M.'^ {d = 2,3) is a polygonal 
domain. Let 7/^ is a quasi-uniform triangulation or rectangular partition of f2 C M'^ with 
mesh size h such that Q = Ui^eTh ^^^o, let (V/^, M^) be a stable mixed finite element 
pair, that is, \h C H^(i7) and C L^(fi) satisfy the inf-sup condition 



(3.1) 



sup 
v,,6Vhnx 



L2 



>/9o||'^h||L2 W^heMhnLliQ), po>0. 



A well-known example that satisfies (3.1) is the following so-called Taylor-Hood ele- 
ment (cf. ^): 



W e C°(fi); v;,|;^ G P2{K) WK e %}, 



In the sequel, we shall only present the analysis for the Taylor-Hood element, but remark 
that the analysis can be readily extended to other stable combinations. On the other 
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hand, constant pressure space is not recommended because that would result in no rate 



of convergence for the approximation of the pressure p (see Section 3.3). 

Approximation space Wh for q variable can be chosen independently, any piecewise 
polynomial space is acceptable as long as Wh D Mh. Especially, Wh C can be 

chosen as a fully discontinuous piecewise polynomial space, although it is more convenient 
to choose Wh to be a continuous (resp. discontinuous) space if Mh is a continuous (resp. 
discontinuous) space. The most convenient and efficient choice is Wh = Mh, which will 
be adopted in the remaining of this paper. 



We now ready to state our semi-discrete finite element method for problem (2.5)-(2.9). 
Letgo = divuo. We seek (u,,,^^, g^) : (0,T] ^ Yhy<MhxWh smd {uh{0),qh{0)) G VfexVr/, 
such that for all t G (0, T] there hold 



(3.2) 


P{Vuh,V\h) - {ph,div\h) 


— (f 7 V/i)dual 




(3.3) 


{divuh,fh) 


= {(lh,^h) 


G Mh, 


(3.4) 


{VUh{0), VWh) = (Vuo, Vw^) , (U;,(0), u) 


= (U07 l^) 


ywh G Vfc, 


(3.5) 




= 


y^h G Wh, 


(3.6) 


{qh{o),Xh) 


= {(lo,Xh) 


VXh G Wh. 



Where qht denotes the time derivative of qh- 

Remark 3.1. (a) We note that ([3^ and (|3^ defines Uh{0) = TZ^Uq and g/i(0) 



Q^qo (see their definitions below). It is easy to see that in general qh{0) 7^ divu^(O) 
although g(0) = divu(O). But it is not hard to enforce qh{0) = divu/j(0) by defining 



u/j(0) slightly differently (in the case of continuous Mh) or simply substituting (3.6) by 



qh{0) '■= divu/i(0) (in the case of discontinuous Mh), even such a modification is not 



necessary for the sake of convergence (see Section 3.3). We also note that Uh{0) = TZ'^Uq 
can be replaced by the projection Uh{0) = Q^uq, the only "drawback" of using the 
projection is that it produces a larger error constant. 

(b) In the case that Mh and/or Wh are discontinuous spaces, since Mh </i H^{Q) and/or 



Wh </i H^{Q), then the second term on the left-hand side of (3.5) must be modified into 
a sum over all elements of the integrals defined on each element K G Th, and additional 
jump terms on element edges may also need to be introduced to ensure convergence. 

We conclude this subsection by citing a few well-known facts about the finite element 
functions. First, we recall the following inverse inequality for polynomial functions [6]: 

(3.7) \\y<Ph\\L^{K)<Coh-^Uh\\LHK) y^hePriK), K eTh. 
Second, for any G L'^{^), we define its projection Q'^cp G Wh as follows: 

{Q''(j),^h) = {<!>, ^h) ^hEWh. 

It is well-known that the projection operator : L^{Q) — Wh satisfies (cf. [5), for any 
G H'{Q) {s > 1) 

(3.8) II QV - <P\\l- + /i||V(QV - 0)|U2 < CH'UWh^, i = min{2, s}. 

We remark that in the case Wh ^ H^{Q), the second term on the left-hand side of the 
above inequality has to be replaced by the broken if^-norm. 

Finally, for any v G H^(i7) we define its elliptic projection T^.'^v G V/^ by 

(V7^V, VW;,) = (VV, VW^) W;, G V;,, 

(7^Vl) = (v,l). 
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Also, for any G H^{Q), we define its elliptic projection S^cp G Wh by 

(V5^'0, V^,) = (V0, V^h) e Wh, 

(5V,1) = (0,1). 

It is well-known that the projection operators TZ^ : H^(i7) — > V/^ and S'^ : H^{Q) Wh 
satisfy (cf. gl E]), for any v G H'{n),(f) G H'{n){s > 1), m = min{3, s} and £ = 
min{2, s} 

(3.9) h~^\\n''^f - ^r\\H-l + ||7^'^v - v||i2 + /^||V(7^'*v - v)|U2 < C/i"'||v||jy,n. 

(3.10) /i-^cSV - 4>\\h-^ + II^V - 4>\\l^ + /ill V(5V - 0)|U2 < Ch'UWH^. 



3.2. Stability and solvability of (|3.2|-([3J1). In this subsection, we shall prove that the 



semi-discrete solution {uh,Ph, qh) defined in the previous subsection satisfies an energy law 



similar to (2.13). In addition, they satisfy the same conservation laws as those enjoyed 



by their continuous counterparts (u,p, g) (see Lemma 2.4). An immediate consequence 
of the stability and the conservation laws is the well-posedness of problem (3.2)-(3.6). 



Moreover, the stability also serves as a step stone for us to establish convergence results 
in the next subsection. 



Lemma 3.1. Let {uh,Ph,qh) be a solution of (3.2)-(3.6) and set ph := Ph + ctQh- Then 
there holds the following energy law: 



(3.11) 

where 



Jhit) + K [ \\ Vph{s) Ilia ds = A(0) Vt G (0,T], 
Jo 



Jh{t) :-- 



Vuh{t) + all qh{t) - 2(f , u/,(t))duai 



Proof. If Uht G L^((0, T); H^(fi)), then (3.11) follows immediately from setting = 
Uht in (3.2), iph = Ph after differentiating (3.3) with respect to t, iph = Ph, and adding the 



resulted equations. On the other hand, if Uht ^ L^((0,T); H^{fl)), then v/j = um is not a 
valid test function. This technical difficulty can be easily overcome by smoothing Uh in t 



through a symmetric mollifier as described in the proof of Lemma 2.3 □ 



Lemma 3.2. Every solution {uh,Ph,qh) of (3.2)-(3.6) satisfies the following conser- 
vation laws: 



(3.12) 



{qh, l) = Cg, {uh, u) = Cu, {ph, l) 



{ph, l) = Cp 



Proof (3.12 )i follows from setting ■i/'/j = 1 in (3.5) and x/i = 1 is (3.6). (3.12)2 follows 
from taking = 1 in (3.3). (3.12)3 derives from letting v/i = x in (3.2). Finally, (3.12)4 
follows from (3.12 )i, (3.12)3, and the definition of ph- □ 

With the help of the above two lemmas, we can prove the solvability of problem 



(13J-(13J. 

Theorem 3.3. Suppose (uo,f) G H\Q) x H-^dQ), and (f , l)duai = 0. then ([3^ 



(3.6) has a unique solution. 

Proof. Since (3.2)-(3.6) can be written as an initial value problem for a system of 
linear ODEs, the existence of solutions follows immediately from the linear ODE theory. 
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To show the uniqueness, it suffices to prove that the problem only has the trivial 
solution if f = uo = 0. For the zero sources, the energy law and the if^-norm stability of 
the elliptic projection. 



Cp 



as 



II Vu;,(0) IU2 = II V7^'^u(0) IU2 < II Vu(0) IU2, 

immediately imply that u,p, and q are constant functions. Since — Cq 
f = Uo = 0, hence, u,p, and q must identically equal zero. The proof is complete. □ 

By exploiting the linearity of equations (3.2)-(3.6), we can prove certain energy esti- 
mates for the time derivatives of the solution {vLh^Vh, qh)- 

Theorem 3.4. The solution {uh,ph,qh) of the semi- discrete method and ph ■= Ph + 
aqh satisfy the following estimates: 

(3.13) Uqh)t\\LHH-r)<Jh{Oy^. 

(3.14) a//3|| Vi{Vuh)t ||l2{l2) + v^ll Vi{qh)t Wl^l^) + v^ll ViVph < 4(0)i 

(3.15) \^\\t{VUh)t\\L°°(L^) + \^\\t{qh)t\\L^{L2) + \^\\t{Vph)t\\L^L^) < [2Jft(0)]i 

(3.16) \\t{qh)tt\\LHH-r) < [2J/.(0)]5. 



Because the proofs of (3.13)-(3.16) follow exactly the same lines of the proofs of their 
differential counterparts given in Theorem 2^ we omit them. 

3.3. Convergence analysis. Define the error functions 

E,(t) := u{t) - Uh{t), Ep{t) := p{t) - phif), Eg{t) := q{t) - q^it), 

and 

p{t) := p{t) + aq{t), ph{t) := ph{t) + aqh{t), Ep{t) := p{t) - ph{t). 
Trivially, Ep{t) = Ep{t) + aEg{t). 



Subtracting each of (3.2)-(3.6) from their respective counterparts in (2.5)-(2.9) we 



obtain the following error equations: 
(3.17) 



(3.18) 
(3.19) 
(3.20) 
(3.21) 



;3(VE„Vv^) - (Ep,divv^) =0 

(divEu,v3/i) = {Eg,(fh) 
u,(0)=7^M0), 
{{Eg)t,^f,)^^^^ + K{\/{Ep + aEg), ViJH) = 

QhiO) = Q'q{0). 



Vv, G V,, 



Introduce the following decomposition of error functions 





= Au + ©u, 


Au 


:= u - 


-7^S 


©u 




- U/x, 


Ep 


= Ap + Qp, 


Ap 


= p- 




Op 


= qV- 




E, 




A. 


= 9 - 


Q\ 




= Q\- 




Ep 


= Ap + Qp, 


Ap 


= p- 


QV, 


ep 


= Q'^p- 


Ph, 


Ep 




% 


■.= p- 


5V, 


% 


= S^P- 


Ph, 



Trivially, <l>p = Ap - ^p + 0p, 6p = 6p + aO^. 
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By the definition of IZ^ and Q!^, tlie above error equations can be written as 

(3.22) ;3(V0u, VV;,) - (0p,divV^) = (Ap,divV;,) Vv^ G V;,, 

(3.23) (div 0u, '^h) = (e„ i^h) - (div A^, y^h) ^^h e M^, 

(3.24) 0^(0) = 

(3.25) ((e,)i, ^h)^,,,, + /t(V$p, V^,) =0 VV'/^ e Wh, 

(3.26) 9^(0) = 0. 

So ©u, 6p, 6p and 6g satisfy tlie same type of equations as Uh,ph,Ph and qh do, except tlie 
terms containing Au and Ap on the right-hand sides of the equations. Hence, ©u, 6p, Qp 
and Qq are expected to satisfy an equahty similar to the energy law (3.11). 

To the end, taking v/^ = {@u)t,Vh = ©p (after differentiating (3.23) with respect to 
t) , iph = = Ap — p + Qp + aQq in (3.22 )-(3.26), and adding the resulted equations we 
obtain 

1 d 



(3.27) 



2dt 



[P\\V&u\\h+a\\Qq\\l.]+4V%\\h 

= {Ap, div (©,)t) - (div (A,)t, Bp) + ((e,)t, ^p - Ap) 

= ^ (Ap, div ©u) + (Og, ^p - Ap) - ((Ap)t, div ©u) 

- (e„ (^p - Ap),) - (div (A Ji, Gp) . 
Integrating in t and using Schwarz inequality we get 

^[/3||V©.(t)||i. + «||0,(t)||y +«:^*||V$p(.)||i.rf. 



(Ap(t),div©,(t)) + (0,(t),^p(t)-Ap(t))- / {eq,{^p{s)-Ap{s))t)ds 

Jo 



((Ap(s))„div©u(s)) + (div(Au(s))„0p(t)) 



ds 



<fliv©u(t) lli. + ^l|Ap(^)lli^ + ^ll0. 



1 



+ -[ii*pWiii^ + iiAp(t)iiy + 2 



2 

L2 



|$p(s)||i.(i5 



+ A/311 V©.(.) Hi. + all 0,(.)||yrf. + - All(^p(^)-Ap(.)).||i.]rf. 
Jo " Jo 



+ C 



(Ap(^)).lli. + («+— )l|div(A„(.)),||i. 



(is. 



Where we have used the facts that ©u(0) = 0, 0g(O) = 0, and 0p = 0p — aQq = 
^p — Ap + <l>p — aQq to get the second inequality, eo is a positive constant to be chosen 
later. Since ($p, 1) = 0, using Poincae's inequality 

II $plU2 < cill V$p ||l2, 
and choosing e = c^^ the above inequality can be written as 

(3.28) /3||V©u(t)||i. + «||0.(t)||i. + K f ||V$p(5)||i. ds 



< / [/3|| V©u(t) 11^2+ all 0,(t)||ycis + ^(t;u,p,p,g), 
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Ait;n,p,p,q) := ||| Ap(t) + ^[\\ ^.(t) + || A,(t) 



+ c 



+ 



a 



^11 (Ap(s)). Hi. + (« + ^)|| div(A,(s)), Hi. 
All (^p(^))«lli^ + ll iMs))s\\l.] ds. 

Jo 



ds 



(3.29) 
(3.30) 

where 



From (3.8)-(3.10) we deduce 

A{t;u,p,p,q) < Ci{T;u,p,p,q) h? , 
A{t;u,p,p,q) <C2{T;u,p,p,q) h^, 



Ci(T;u,p,p,g) ■.= C[/3-'\\p\ 



+ P ^Wvt IIl2(//i) 



+ (a + ciK ^)|| (divu)t 11^2(^1) + a ^ [11 :P lli-(i?i) + lb* IIl2{hi)] ^ 



C2(T;u,p,p,g) :=C[rip| 



2 , II ~ ||2 

L°°(H2) + P II I't IIl2(//2) 



+ (a + ciK ^) II (divu)t 11^2(^2) + a ^\\\p\\1^(^h^) + WPtWh^n^)]. 



It follows from an application of the Gronwall's lemma to (3.28) that 

(3.31) max[/?||Veu(t)||i2 + a||e,(t)||y +«: / \\V%{s)\\l. ds < A{T-xy,p,q) . 

— — •/ 

An application of the triangle inequality yields the following main theorem of this 
section. 

Theorem 3.5. Assume uq G H^(r2). Let 

dl{T; U,p,q) := C [/?|| U ||ioo(j:^2) + a|| g ||ico(^l) + k|| P 11^2(^:^2)] , 

C2{T; u,g) := C[/?|| u ||ioo(j:^3) + a|| g ||ico(^2)] . 
Suppose that Ci{T; u,p,p, q) < oo and Ci{T; u,p, q) < oo, then there holds error estimate 

(3.32) max^[V^|| V(u(t) - u,(t))|U2 + v^||g(t) - g,(t)|U2] 

+ [k \\V{p~pH)mhdtY < [Ci(T;u,p,p,g)ie?+5i(T;u,p,g)^]/i, 

Moreover, if C2{T;u,p,p,q) < oo and C2{T;u, q) < oo, then there also holds 

(3.33) max [v^||V(u(t) - u,(t))|U2 + Mm - QhimL^] 

< [C2iT; u,p,p, g)i e? + C2{T; u, q)^ h\ 



Remark 3.2. (a) We note that the above error estimates are optimal, and || V$q ||l2(l2) 
II V(p - S^p) II ^2(2^2) enjoys a superconvergence when the PDE solution is regular. It is 
also interesting to remark that the initial errors at t = do not appear in the above error 
bounds. 

(b) In light of Theorem 2.6, the regularity assumptions of Theorem 3.5 are valid if the 
domain Q and datum functions f and Uq are sufficient regular. 
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4. Fully discrete finite element methods. 

4.1. Formulation of fully discrete finite element methods. In this section, we consider 
space-time discretization which combines the time-stepping scheme of Algorithm 1 with 
the (multiphysical) spatial discretization developed in the previous section. We first prove 
that under the mesh constraint At = OQi^) the fully discrete solution satisfies a discrete 
energy law which mimics the differential energy laws (2.13) and (2.11). We then de- 



rive optimal order error estimates in various norms for the numerical solution, which, as 
expected, are of the first order in At. 



Using the time-stepping scheme of Algorithm 1 to discretize (3.2)-(3.6) we get the 



following fully discrete finite element method for problem (2.5)-(2.9) 
Fully discrete version of Algorithm 1: 

(i) Compute g° G Wh and u° G V/^ by 

(4.1) {ql.Xh) = {qo,Xh) 

(4.2) (Vu° , Vw;,) = ( Vuo, Vw^) , (u° , v) = (uq, v) 

(ii) For = 0, 1, 2, ■ ■ ■ , do the following two steps 
Step 1: Solve for eV^xMh such that 



Vw/, G V,,, 



(4.3) 
(4.4) 
(4.5) 



/3(V<+\VV,) - (p;:+\divV,) = (f,V,)dual 



IV u 



Jl+1 
h ■ 



^h) = {qh, Vh) 



Step 2: Solve for g^+^ G Wh such that 

(4.6) (rftgrS^.) + «:(V(agr^ -Fp^+^), VtA.) = 

(4.7) (gr\l)=C,. 



V^, G Wh, 



Some remarks need to be given before analyzing the algorithm. 



Remark 4.1. (a) At each time step, problem (4.3)-(4.5) in Step 1 of (ii) solves a 
generalized Stokes problem with a slip boundary condition for u. Two conditions in (4.5) 



are imposed to ensure the uniqueness of the solution. Well-posedness of the Stokes problem 



follows easily with help of the inf-sup condition (3.1). 

(b) g^^"^ is clearly well-defined in Step 2 of (ii). 

(c) The algorithm does not produce (and p\), which is not needed to execute the 
algorithm. 

(d) Reversing the order o/ Step 1 and Step 2 in (ii) of the above algorithm yields the 
following alternative algorithm. 

Fully Discrete Finite Element Algorithm 2: 

(i) Compute g° G Wh and (u°,]5^) G V/^ x M/^ by 

{qt Xh) = (go, Xh) ^Xh e Wh. 

f3{Vul VV,) - (pO, div V,,) = (f , V,,)dual Vv, G V;,, 

(divu^, ifh) = {qh, ^h) e Mh, 

(pO, 1) = Cp, (u°, v) = C^. 
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(ii) For n = 0,1,2, ■■■ , do the following two steps 
Step 1: Solve for q^^^ G Wh such that 

Step 2: Solve for {yLl+\^+^) eVhxMh such that 



/3(Vu;:+\ VV,) - (p;j+\divV,) = (f,V,)dual 



V</^/^ G Mh, 



1) 



We note that Algorithm 2 requires the starting values q^ and {vl\,jP^, and the latter is 
generated by solving a generalized Stokes problem att = 0. In general, u° of Algorithms 
1 and 2 are different. Later we will remark that Algorithm 2 also provides a convergent 
scheme. 

4.2. Stability analysis of fully discrete Algorithm 1. The primary goal of this subsection 
is to derive a discrete energy law which mimics the PDE energy law (2.11 ) and the semi- 
discrete energy law (3.11). It turns out that such a discrete energy law only holds if h 
and At satisfy the mesh constraint At = 0{h'^). This mesh constraint is the cost of using 
the time delay decoupling strategy in (4.4). 

Before discussing the stability of the fully discrete Algorithm 1, We first show that 
the constraints (4.5) an ( |4.7 ) are consistent with the equations (4.3), (4.4), and (4.6). 

Lemma 4.1. Let {{vLl,pl,ql)}n>i satisfy (|43|, (|0|), and iQ. Define pi := p^ + 
aq]l~^ for n = 1,2,3, ■■■ . Then there hold 



(4.8) 
(4.9) 
(4.10) 
(4.11) 



Cp 



for n 
for n 
for n 
for n 



0,1,2, 
1,2,3, 
1,2,3, 
1,2,3, 



where Cg,Cu,Cp, and Cp are defined by (1.14)-(1.16) and (1.18) 
Proof. Taking i/'/i = 1 in (4.6) yields 



(«+\l)=0. 



Hence, 



l) = (g°,l) = (go,l)=C, 



for n = 0,1,2, 



So (4.8) holds. (4.9) immediately follows from setting (ph = 1 in (4.4) and using (4.8). 

As in the time-continuous case, (4.10) follows from taking v/j = x in (4.3) and using 
the fact that div x = d and Vx = /. 

Finally, (4.11) is an immediate consequence of the definition p^ := p^ + aq^^ , (4.8), 
and (4.10). The proof is complete. □ 

Next lemma establishes an identity which mimics the continuous energy law for the 
fully discrete solution of Algorithm 1. 

Lemma 4.2. Let {(u^,p^, gj^)}„>i be a solution of ([4^1]) -([ij)) . Define p^ := pl+aq]^'^ 
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for n > 1. Then there holds the following identity: 

(4.12) j-r + ^tj2[4 vj>r^ wh + dt'^^' wi^ + "II « iiio 



n=0 



n=0 



where 
(4.13) 



7^ - 



■||v<+Mli. + «lk^lli^-2(f,<+i)d.ai 

Proof Setting v/j = dt'n^'^'^ in (4.3) gives 



0,1,2,- 



(4.14) 



Applying the difference operator dt to (4.4) and followed by taking the test function 
= K^' = Pi"-' - aql yield 

(4.15) (divd,ur\p;:+i) = {d,ql,pl^') - a{dtqtql) 

a 



{dtql,pt')-^[dt\\ql\\h+mdtql\\l. 



To get an expression for the first term on the right-hand side of (4.15 ), we set ip^ = p^ 



71+1 



in (4.6) (after lowering the index from n + 1 to n in the equation) and using the definition 
(4.16) 



pl'-^K^' + aq^ito get 



Vpl-^'Wl + nAtid.VK^Wpl-''). 



-K\ 



(4.17) 



Combining (4.14)-(4.16) we obtain 
1 



-d^ 



Vu^i Hi. + all g," Hi. - 2(f , nl-^\^^,\ + «:|| Vp^^ ||i. 
At 



+ 



^tVu^^ Hi. + all d,q]: Hi. = -^At{d,Vp^,+\Vpr')). 



Applying the summation operator At Yln=o ^ — ^ — 'Kt ~ ^ yields 



(4.18) 



^n+l II 2 

Il2 



At 



dtVnl-^' lli. + «ll«^' IliO 



n=0 



- ^{AtrJ^idtWt^'^Pr'). 



n=0 



where is defined by (4.13). Hence, (4.12) holds. The proof is completed. □ 

Remark 4.2. The extra term on the right-hand side of (4.12) is due to the time lag 

in (4.4) which is needed to decouple the original whole system into two sub-systems. This 

term is smaller if high order time-stepping schemes are used to replace the implicit Euler 

scheme in Algorithm 1. 

Theorem 4.3. Let {(uJJ,p^, g)^)}n>i he same as in Lemma 4-2. Set q^^ := Then 

there holds the following discrete energy inequality: 



(4.19) jr+^tY.l'^wvp- 

n=0 
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provided that At = 0{h'^). 

Proof. By Schwarz inequality and the inverse inequality we have 



(4.20) \{dNrK^\Vpt')\ 



< 



At 
1 

At 1 



Il2 + 



h 1IL2 
n+1 ||2 



r>2/)-2|| ^i^+l _ ||2 _|_ II Vti' 
Cq"- WPh + -^^ "^^ 



To bound the first term on the right-hand side, we appeal to the inf-sup condition (3.1) 
(note that p^+^ ~pleLl{VL)) and ^ to get 

(divv,,p):+^-p^) 



(4.21) 



sup 



Po ^ sup 



II Vv/, 11^2 
/3(V(uri-0,Vv,) 



Vv 



h ||l2 



1^2 



/3o-VAtN*V<+HU2. 



( |4.19[ ) follows from combining ( |4.12[ ), ( |4.20[ ), and ( |4.21[ ) and choosing At < -^^fi^- 
The proof is complete. □ 

Remark 4.3. It can be shown that all the results obtained in this subsection for 
Algorithm 1 also hold for Algorithm 2. The main differences are (a) the "correct" discrete 
pressure p^ for Algorithm 2 is p^ := p^ + aq"^ instead ofp^ := Ph + cuQh'^ ! (b) the "correct" 
energy functional J'^ for Algorithm 2 is 

^ 11^2 + all glllia -2(f,4)dual 



0,1,2, 



4.3. Convergence analysis. In this subsection we derive error estimates for the fully 
discrete Algorithm 1. To the end, we only need to derive the time discretization errors 
by comparing the solution of (4.1)-(4.7) with that of (3.2)-(3.6) because the spatial 
discretization errors have already been derived in the previous section. 

Introduce notations 



:= u,(t„) - ul 

Ph{tn) 



It is easy to check that 



E; = E; + aAtdtqhitn) 



e: 



Ph{tn) 



vl- 



El + aE^ 



n-\ 



Lemma 4.4. Let {(u^,p)J, g}^)}„>o be generated by the fully discrete Algorithm 1 and 



E'^,Eg,E~, and E^ be defined as above. Set qh{t- 
holds the following identity: 



:= qh{to) and := Then there 



(4.22) 



+ At 



n 

n=0 



ve;+' Hi. + 



At 



dtVY.l+^\\l2 + a\\dtE: 



Ii2) 



- (At)2^[«:(rf,VE|-^\ V^r^) - (ci,V(Wi),i?r') 

ra=0 

+ AtY,{RtE;^'), 



n=0 



FINITE ELEMENT APPROXIMATIONS OF SWELLING DYNAMICS OF GELS 



19 



where 

(4.23) 
(4.24) 



1 r 
2 



vypi^+l ||2 I ^11 Tp(. ||2 
V-t'u IIl2 + "II IIl2 



^h-=-^. {s-t,^i){qh)tt{s)ds. 

J tn — l 



0,1,2, 



Proof. First, by Taylor's formula and (3.5) we have 



(4.25) {dtqh{tn+i),A) + 4V{ph{tn+i) + aqh{tn+i)),V^h) = {RT\^h) "^^h e Wh, 



Then, subtracting (4.3) from (3.2), (4.4) from (3.3), (4.5) from (3.12), and (4.6) from 



(4.25), (4.2) from (3.4) respectively, we get the following error equations: 



(4.26) 
(4.27) 
(4.28) 
(4.29) 
(4.30) 



/5(VEr\ Vv,) - (i?|+\divv,) = 
{EZ^\u) = 0, (E|+\1)=0, E°=0, 

(e;+\i) = o, eJ = o. 



( [422] ) then follows from setting v,, = in ( [426] ), (ph = E~^^ in ( [4271 ) (after 



applying the difference operator dt to the equation), iph = Ep~^^ = E~^ + aE"' in (4.29) 



(after lowering the index from n + 1 to tt, in the equation), adding the resulted equations 
and applying the summation operator At Yln=o sum. The proof is complete. □ 

From the above lemma we then obtain the following error estimate. 
Theorem 4.5. Let N be a (large) positive integer and At :— ^. 



Let the error 



functions E'^,E^,E~, and i?" be same as in Lemma 4-4 Assume At satisfies the mesh 



O 9 

constraint At < „ ° 2 ^ • Then there holds error estimate 



N-1 



(431) 



max 

o< 



lax^fv^ll VE;^ IU. + v/^ll E:^ IU.] + (At J2 '^11 V^;^' 

n=0 

N-1 1 

[E (^11 ^(Eu^' - K) + «ll - E^' IliO] ' < CsiT; g.)^ At, 



n=l 



where 



CiT; g,) := 16[/3^'P\\ (g,)t ||i.(^.) + k'^ (q^h 



L2(/f-l) 



]■ 



Proof. To derive the desired error bound, we need to bound each term on the right- 
hand side of (4.22). Before doing that, on noting that E'^ = and E° = and (4.23) we 
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rewrite (4.22) as 



(4.32) £t' + AtK\\VEl\\l, + ^\\El\\l 



n=l 



At 



+ AtJ2M VE;^' Hi. + — (/3|| rf^VEr^ Hi. + all d,E^ ||i.) 



-{AtrJ2[<dtVE^+\VE;^') - Kg.(Wi),i?f ^) 

n=0 

e 

+ AtJ2{RlE;^')- 



n=0 



We now estimate each term on the right-hand side of (4.32). First, by Schwarz in 



equahty, the inverse inequahty (3.7), the inf-sup condition (3.1), and (4.26) we have 



(4.33) k\ {dtVE^+\ VE;+i) I < Kcoh-^ VE;+^ \\lA\ dtEI+' 

(divVfe,ci,E|+^) 

livv.iu. 



< KCoh~^(3oA\VE^+'\\L2 sup 



<KCoh-'Po \\'^E;^'\\l2 sup ,, 

v^eVfe ||Vv/i||l2 

< KCoh~'(3^'P\\ VE;+i IU^II dtVEl+' 



To bound the second term on the right-hand side of (4.32), we first use the summation 
by parts formula and dtquito) = to get 

(4.34) 5^(rf2g,(t„_,,),E^+i) = Uqu{t,^,),E'^') - Y,{dMtn) , d,E^+') . 



n=0 



n=l 



We then bound each term on the right-hand side of (4.34) as follows: 

II VTT7i^+l l|2 I o-2n\\ \ ||2 



(4.35) 



< 



4(^^)2 II ^^u^ 11^2 + ^0 /5|l i(lh)t \\L2{{tt,ti+iy,L^)- 



(4.36) 



n=l 



J2{dtqhitn),dtE^+') \ < II dtQkitn) \\l4 dtE^+^ 1^. 

n=l 

< f3o'(3j2\\ dtqhitn) \\l4 rf^VE^' ||l2 



ra=l 



< 



Ell^*^ErMli^+/?o-'/?ll i<lH)t 



Il2(L2)- 



n=l 
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Finally, we bound the third term on the right-hand side of (4.32) by 

(4.37) \{RIe;^')\ < \\K\\H-^\\VE;+'h^ 

< 



-\\ V^p+^ 11^2 + — II {qh)tt Ili2({t„_i,t„);ff-1), 



where we have used the fact that 



^„||2 



'^h \\H 



-1 < 



At 



{<lh)tt{t) ll^-i dt. 



Substituting (4.33)-(4.37) into (4.32) and using the assumption that At < 



yields 



E: 



-1 ||2 



l.+a\\El\\l.+M[.\\VEl\\l. + '^\\E]\\l.] 

+ At^[«:|| V^;+i lli^ + — (/?|| (i,VEr^ lli^ + all d,El ||iO 
< 16(At)2[/3o-'/?|| (g.)* 11^2(^2) + iqnh Wl^^H-r)], 



n=l 



which trivially implies (4.22). The proof is complete. □ 
Th 

error t 

(4.38) 



Theorem 4.6. The solution of the fully discrete Algorithm 1 satisfies the following 
error estimates: 



max Jv//3|| V(u(t„) - <) ||i2 + v^|| g(t„) - ql H^a] 



0<n<Af 



N 



At5^«:||V(p(tn)-K) 



Il2 



n.=0 



< [Ci(T; u,p,p, g)i e? + C,{T- u,p, g)i] /i + C;{T- q^^At, 

provided that Ci(T; u,p,p, q) < oo, Ci(T; u,p, q) < oo, Cs(T; qA < oo, and At < „ ^° 2 • 

Moreover, if, in addition, C2{T;u,p,p, q) < 00 and C2{T;u, q) < 00, then there also 
holds 



(4.39) 



max [v^ll V(u(t„) - ul) \\l2 + v^|| g(t„) - g^J 1^2] 

0<n<A'' 

< [C2(T; u,p,p, g)i e? + C^iT; u,q)^h^ + C^{T; q^f^At. 
Proof. The assertions follow easily from first using the triangle inequality on 

U{tn) - < = Eu{tn) + El, 

q{t^)-ql = E,{t^) + E-^, 
p{tn)-pl = E,{tr,) + E;, 

and then appealing to Theorems 4.5| and 4^ □ 

Remark 4.4. (a) In light of Theorems 2.6 and 8.4, the regularity assumptions of 
Theorem \4.()\ are valid if the domain Q and datum functions f and uq are sufficient regular. 

(b) It can be shown that all the results proved in this subsection for Algorithm 1 still 
hold for Algorithm 2. The main differences are (i) the "correct" discrete pressure pf^ for 
Algorithm 2 is p'^ '■= Ph ~^ '^^h> (V "correct" error functional for Algorithm 2 is 

1 



VE^J|i2+a||E ||i2 



0,1,2,--- 
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5. Numerical experiments. In this section we present some 2-D numerical experiments 
to gauge the efficiency of the fully discrete finite element methods developed in this paper. 
Three tests are performed on two different geometries. The gel used in all three tests is 
the Ploy(N-isopropylacrylamide) (PNIPA) hydrogel (cf. [16] and the references therein). 
The material constants/parameters, which were reported in [16], are given as follows: 



E = Qx 10^ 
u = 0.43, 

K " 



G 



3(1 -2z/) 
E 

2(1 + v) 



14285.7, 



2097.9, 



Young's modulus, 
Poisson's ratio, 

bulk modulus, 
shear modulus. 



Two other material constants/parameters, which were not given in [16], are taken as 
follows in our numerical tests: 

(p = 0.15, porosity, 

^ = 100, friction constant. 

In addition, we use the following initial condition in all our numerical tests: 

Uo(x) = 10"'^ sin(a;i + X2)(l, 1). 



Test 1: Let Q = (0, 1)^. The external force is taken as 

f =(/l,/2)=0.1t„, 

where tt„,g„„t denotes the unit (clockwise) tangential vector on dQ. Note that the compat- 
ibility condition 



(5.1) 

is trivially fulffiled. 



/ i{x)ds = o 

Jan 



Fig. 5.1. Computational domain, boundary data, and mesh of Test 1 



Figure 5J^ shows the computational domain, color plot of the force function /i + /2, 
and the mesh on which the numerical solution is computed. The mesh consists of 2360 
elements, and total number of degrees of freedom for the test is 12164. At = 0.01 is used 
in this test. 



Figure |5.2| displays snapshots of the computed solution at three time incidences. Each 
graph contains color plot of the computed pressure '■= Ph + c^Qh'^ and arrow plot 
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of the computed displacement field u]^. The three graphs on the first row are plotted 
on the computational domain Q = (0, 1)^, while three graphs on the second row are 
respectively deformed shape plots of the three graphs on the first row with 500 times 
magnification, which shows the deformation of the square gel under the mechanical force 
f on the boundary. As expected, the gel is slightly rotated clockwise and is slightly bent 
near the top and bottom edges. We also note that the expected conserved quantities are 
indeed conserved in the computation, their respective values are given as follows: 

Cu= dS = 9.935 X 10"^ Cg = q]l{x) dx = 9.935 x 10~^ 

Jan Jq 



Cp= I pI{x) dx = 1.489. 
Jn 






Fig. 5.2. Test 1: color plots of computed pressure pJJ and arrow plots of computed displacement uJJ ai n = 0, 3, 5 (first 
row). Deformed shape plots of the graphs on the first row (second row). At = 0.01. 



Test 2: This test is same as Test 1 except f = (/i, is changed to 

0.5 for xi = 0, < X2 < 1, 
/^(a;) := <j -0.5 for Xi = 1, < < 1, 

for X2 = 0, 1, < xi < 1, 



/2 



X 



:= 0. 



So parallel forces of opposite directions are applied at the left and the right boundary of 



the square gel. Clearly, the compatibility condition (5.1) is satisfied 



Like Figure 5J^, Figure 5^ shows the computational domain, color plot of the force 
function /i + /a- The mesh parameters are also same as those of Test 1, including the 
time step At = 0.01. 



Figure 5^ is the counterpart of Figure 5^ for Test 2. Since parallel forces of equal 
magnitude but opposite directions are applied on the left and the right boundary of the 
square gel, so the gel is squeezed horizontally. Because of the incompressibility of gel, the 
total volume must be conserved. As a result, the deformation in the vertical direction 
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Fig. 5.3. Computational domain, boundary data, and mesh of Test 2 





I I i ^ I ^ 
111/ J / 





FlG. 5.4. Test 2.- co/or plots of computed pressure p^ and arrow plots of computed displacement uJJ ai n = 0, 3, 5 (first 
row). Deformed shape plots of the graphs on the first row (second row). At = 0.01. 



is expected. The second row of Figure 5.4 precisely shows such a deformation (with 500 



times magnification). It should be noted that the swelling dynamics of the gel goes super 
fast, it reaches the equilibrium in very short time. For the conserved quantities, Cq and 
Cu are same as those in Test 1, and the other two numbers are given by 



Cp= [ p^{x) dx = 0.321, Cp= [ pI{x) dx = 1.809. 
Jq Jn 



Recall that Cp and Cp depend on the force function f . 

Test 3: Same material parameters/constants and initial condition uq as in Tests 1 
and 2 are assumed. However, the computational domain is changed to the following one 



Q:=\xeW 



X 1 



+ 



0.16 0.04 - " 



and the external force function f is taken as 



0.5 for X e dn, -0.2 < xi < 0, 
0.5 for X e dQ, < xi < 0.2, 



Mx) - 
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As in Test 2, the above f means that parallel forces of same magnitude but opposite 
directions are applied at the left half and the right half ellipse (boundary), however, these 
two forces now collide at two points (0, ±0.2) on the boundary. Clearly, the compatibility 
condition (5.1) is satisfied. 






4I> 



<1> 



Fig. 5.6. Test 3: color plots of computed pressure pJJ and arrow plots of computed displacement uJJ at n = 0, 1, 3 (first 
row). Deformed shape plots of the graphs on the first row (second row). Ai = 0.001. 



Figure 5.5 is the counterpart of Figures 5.1 and 5.3 for Test 3. The mesh consists of 
1200 elements, and total number of degrees of freedom for the test is 6244. A smaller 
time step At = 0.001 is used in this test. 

Figure |5.6 is the counterpart of Figures 5.2 and 5.4 for Test 3. The deformations 
displayed on the second row are magnified by 20 times instead of 500 times as in Figures 



5.2 and 5^ Due to curve boundary, the applied parallel forces of same magnitude but 
opposite directions collide at the boundary points (0, ±0.2). It is expected that the gel 
should buckle around those two points, which is clearly seen in the plots on the second row 



of Figure |5.6[ We also note that because relatively bigger forces are applied on the gel, 
the swelling dynamics of the gel goes even faster, hence, reaches the equilibrium quicker. 
This is the main reason to use a smaller time step for the simulation. 
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The conserved quantities for Test 3 are given as follows: 



Cu = / ul{x) dS = 4.902 X 10-^ Cg^ [ ql{x) dx 



4.902 X 10 



-6 




0.178. 
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